#Set seed for reproducibility
set.seed(123)
# Parameters
n_sim <- 10000
D1 <- 2
# Simulate r and g
r <- rnorm(n_sim, mean = 0.08, sd = 0.01)
g <- rnorm(n_sim, mean = 0.03, sd = 0.005)
# Compute values
value <- D1 / (r - g)
# Filter for valid values where r > g
value <- value[is.finite(value) & (r > g)]
# Plot
library(ggplot2)
library(scales)
ggplot(data.frame(Value = value), aes(x = Value)) +
geom_histogram(bins = 60, fill = "#2c7fb8", color = "white", alpha = 0.8) +
scale_x_continuous(labels=scales::dollar)+
labs(
title = "Monte Carlo Valuation Distribution (Gordon Growth Model)",
x = "Firm Value",
y = "Frequency"
) +
theme_minimal(base_size=20)Equity Valuation and Simulation
Monte Carlo Implementation, continued
Exercise 1: changing revenue growth
#Set seed for reproducibility
set.seed(123)
# Initialize an empty dataframe
results <- numeric(0)
sim_assumptions <- netscape_assumptions
n_sim=10000
#Run the simulation
for (i in 1:n_sim){
# Simulate revenue growth for 10 years from a normal distribution
rev_growth <- rnorm(n_years, mean = 0.65, sd = 0.1)
# Modify assumptions with simulated growth
sim_assumptions$revenue_growth <- rev_growth
# Compute EV for this simulation
ev <- get_EV(sim_assumptions,
netscape_base_rev,
sim_assumptions$terminal_r,
sim_assumptions$terminal_growth)
results[i] <- ev
}
#Create a ggplot histogram
results%>%
as.tibble()%>%
ggplot(aes(x = value)) +
geom_histogram(fill = "skyblue", color = "white", bins = 50)+
scale_x_continuous(labels = scales::dollar) +
labs(
title = "Monte Carlo Simulation of Enterprise Value",
subtitle = paste0('Drawing ',comma(n_sim), ' simulations of Revenue Growth.'),
x = "Enterprise Value (USD, Thousands)",
y = "Frequency"
) +
theme_minimal(base_size = 14)
Exercise 2: Changing Revenue Growth and COGS
#Set seed for reproducibility
set.seed(123)
# Initialize an empty dataframe
results <- numeric(0)
sim_assumptions <- netscape_assumptions
n_sim=10000
#Run the simulation
for (i in 1:n_sim){
# Simulate revenue growth for 10 years from a normal distribution
rev_growth <- rnorm(n_years, mean = 0.65, sd = 0.1)
cogs_perc <- rnorm(n_years, mean = 0.10, sd = 0.025)
# Modify assumptions with simulated growth
sim_assumptions$revenue_growth <- rev_growth
sim_assumptions$cost_of_sales_pct <- cogs_perc
# Compute EV for this simulation
ev <- get_EV(sim_assumptions,
netscape_base_rev,
sim_assumptions$terminal_r,
sim_assumptions$terminal_growth)
results[i] <- ev
}
#Create a ggplot histogram
results%>%
as.tibble()%>%
ggplot(aes(x = value)) +
geom_histogram(fill = "skyblue", color = "white", bins = 50)+
scale_x_continuous(labels = scales::dollar) +
labs(
title = "Monte Carlo Simulation of Enterprise Value",
subtitle = paste0('Drawing ',comma(n_sim), ' simulations of Revenue Growth and % COGS.'),
x = "Enterprise Value (USD, Thousands)",
y = "Frequency"
) +
theme_minimal(base_size = 20)
Exercise 3: Changing Revenue Growth, COGS, and CAPEX
#Set seed for reproducibility
set.seed(123)
# Initialize an empty dataframe
results <- numeric(0)
sim_assumptions <- netscape_assumptions
n_sim=10000
#Run the simulation
for (i in 1:n_sim){
# Simulate revenue growth for 10 years from a normal distribution
rev_growth <- rnorm(n_years, mean = 0.65, sd = 0.1)
cogs_perc <- rnorm(n_years, mean = 0.10, sd = 0.025)
capex_perc <- seq(rnorm(1, mean = 0.45, sd = 0.1),to=0.1,length.out = 10)
# Modify assumptions with simulated growth
sim_assumptions$revenue_growth <- rev_growth
sim_assumptions$cost_of_sales_pct <- cogs_perc
sim_assumptions$capex_pct <- capex_perc
# Compute EV for this simulation
ev <- get_EV(sim_assumptions,
netscape_base_rev,
sim_assumptions$terminal_r,
sim_assumptions$terminal_growth)
results[i] <- ev
}
#Create a ggplot histogram
results%>%
as.tibble()%>%
ggplot(aes(x = value)) +
geom_histogram(fill = "skyblue", color = "white", bins = 50)+
scale_x_continuous(labels = scales::dollar) +
labs(
title = "Monte Carlo Simulation of Enterprise Value",
subtitle = paste0('Drawing ',comma(n_sim), ' simulations of Revenue Growth, % COGS, and % CAPEX.'),
x = "Enterprise Value (USD, Thousands)",
y = "Frequency"
) +
theme_minimal(base_size = 20)
final_results=results
Exercise 4: Different levels for terminal \(r\) and \(g\)
#Set seed for reproducibility
set.seed(123)
# Initialize an empty dataframe
results <- numeric(0)
new_data <- data.frame()
sim_assumptions <- netscape_assumptions
n_sim=1000
g_sequence <- seq(0.05,0.03,length.out=10)
r_sequence <- seq(0.15,0.225,length.out=10)
for(s in 1:10){
#Run the simulation
for (i in 1:n_sim){
# Simulate revenue growth for 10 years from a normal distribution
rev_growth <- rnorm(n_years, mean = 0.65, sd = 0.1)
cogs_perc <- rnorm(n_years, mean = 0.10, sd = 0.025)
capex_perc <- seq(rnorm(1, mean = 0.45, sd = 0.1),to=0.1,length.out = 10)
# Modify assumptions with simulated growth
sim_assumptions$revenue_growth <- rev_growth
sim_assumptions$cost_of_sales_pct <- cogs_perc
sim_assumptions$capex_pct <- capex_perc
sim_assumptions$terminal_growth = g_sequence[s]
sim_assumptions$terminal_r = r_sequence[s]
# Compute EV for this simulation
ev <- get_EV(sim_assumptions,
netscape_base_rev,
sim_assumptions$terminal_r,
sim_assumptions$terminal_growth)
results[i] <- ev
}
#Store the Pairs
new_data=new_data%>%
rbind(data.frame(Scenario = paste0('Scenario ', s),
EV=results))
message(paste0('Finished simulation for pair number ',s,'.'))
}
#install.packages('ggridges')
library(ggridges)
#Create a ggplot histogram
new_data%>%
mutate(Scenario = factor(Scenario,levels=paste0('Scenario ',10:1)))%>%
ggplot(aes(x = EV,
y = Scenario,
fill=..x..)) +
geom_density_ridges_gradient(scale = 3, rel_min_height = 0.01)+
scale_x_continuous(labels = scales::dollar) +
scale_fill_viridis_c(name = "Enterprise Value", option = "plasma",labels=scales::dollar)+
labs(
title = "Monte Carlo Simulation of Enterprise Value",
subtitle = paste0('Drawing ',comma(n_sim), ' simulations, best-to-worse scenarios.'),
x = "",
y = "Frequency",
fill = "Enterprise Value (USD, Thousands)")+
theme_minimal(base_size = 14)+
theme(legend.position = 'bottom',
legend.title = element_text(size = 14),
legend.text = element_text(size = 12),
legend.key.height = unit(0.6, "cm"),
legend.key.width = unit(2.5, "cm"))
Final Question: should we invest?
#Create a ggplot histogram
prices=final_results%>%
as.tibble()%>%
mutate(Price=value/netscape_assumptions$shares_outstanding)%>%
mutate(Situation=ifelse(Price>=28.37,'Undervalued','Overvalued'))
prices%>%
ggplot(aes(x = Price)) +
geom_histogram(aes(fill = Situation), bins = 50)+
scale_x_continuous(labels = scales::dollar)+
scale_fill_manual(values=c('darkgreen','darkred'))+
geom_vline(xintercept=28.37,linetype='dashed',size=1)+
annotate(geom='text',x=40,y=500,label= paste0('True Price is higher than \nbaseline price only in ',
percent(mean(prices$Price>=28.37)),
' of the cases!'))+
labs(
title = "Monte Carlo Simulation of Share Price",
subtitle = paste0('Drawing 10,000 simulations'),
x = "Enterprise Value (USD, Thousands)",
y = "Frequency"
) +
theme_minimal(base_size = 20)+
theme(legend.position='bottom')
Base Case Simulation
- The base case denotes the Free Cash Flow in period \(t\) as:
\[ FCF_t = Sales_t\times(1-COGS) \pm CAPEX \pm \Delta NWC \]
- Initial Sales levels are \(\$100\)
- The percent of Cost of Goods Sold (COGS) is \(60\%\)
- Capital Expenditures (CAPEX) is \(15\%\) of Sales
- Net Working Capital investment is \(5%\) of sales
library(glue)
library(MASS)
set.seed(123)
n_sim <- 10000
# Parameters
initial_sales <- 1000
cogs_pct <- 0.6
capex_pct <- 0.15
nwc_pct <- 0.05
simulate_values <- function(r, g, label='All') {
firm_value <- numeric(n_sim)
for (i in 1:n_sim) {
sales <- initial_sales
nwc_last <- sales * nwc_pct
fa_last <- sales * capex_pct
fcf <- numeric(5)
for (t in 1:5) {
sales <- sales * (1 + g[i])
nwc_current <- sales * nwc_pct
fa_current <- sales * capex_pct
delta_nwc <- nwc_current - nwc_last
delta_capex <- fa_current - fa_last
nwc_last <- nwc_current
fa_last <- fa_current
fcf[t] <- sales * (1 - cogs_pct) - delta_capex - delta_nwc
}
fcf_terminal <- fcf[5] * (1 + g[i])
perp_value <- fcf_terminal / (r[i] - g[i])
discounted_fcf <- sum(fcf / (1 + r[i])^(1:5))
discounted_perp <- perp_value / (1 + r[i])^5
firm_value[i] <- discounted_fcf + discounted_perp
}
firm_value <- firm_value[is.finite(firm_value) & (r > g)]
data.frame(Value = firm_value, Scenario = label)
}
# Normal distributions
r_norm <- rnorm(n_sim, mean = 0.08, sd = 0.01)
g_norm <- rnorm(n_sim, mean = 0.03, sd = 0.005)
df_norm <- simulate_values(r_norm, g_norm, "Normal")
base_norm <- simulate_values(r_norm, g_norm, "Normal")
ggplot(base_norm, aes(x = Value)) +
geom_histogram(bins = 60, fill = "#2c7fb8", color = "white", alpha = 0.8)+
scale_x_continuous(labels = dollar) +
labs(
title = "Monte Carlo Valuation Distribution Normal vs. Alternative Distributions",
subtitle = glue('Based on {comma(n_sim)} simulations.'),
x = "Firm Value",
y = "Frequency",
fill = "Scenario",
color = "Scenario"
) +
theme_minimal(base_size = 20)+
theme(legend.position = 'bottom')
Extension 1: Alternative Distributions
Key Point: Monte Carlo doesn’t assume normality — any distribution that fits the problem can be used!
While normal distributions are common, other distributions (e.g., uniform, log-normal, etc.) can also be used:
For example, if \(r \sim \text{Uniform}(0.05, 0.10)\), the growth rate could be chosen randomly from a uniform distribution between \(5\%\) and \(10\%\)
Distributions like the log-normal may be more appropriate for modeling returns, as they respect the non-negativity constraint - $r $. Therefore, for a given variable \(X\), it could be sampled from any distribution
In what follows, we’ll simulate results drawing from an Uniform distribution for \(r\) and a Beta for \(g\)
# Baseline case
r_norm <- rnorm(n_sim, mean = 0.08, sd = 0.01)
g_norm <- rnorm(n_sim, mean = 0.03, sd = 0.005)
df_norm <- simulate_values(r_norm, g_norm, "Normal")
# Alternative distributions
r_unif <- runif(n_sim, min = 0.06, max = 0.10)
g_beta <- 0.06 * rbeta(n_sim, 2, 5)
df_alt <- simulate_values(r_unif, g_beta, "Alternative")
# Combine results
df_all <- rbind(df_norm, df_alt)
# Plot combined chart
ggplot(df_all, aes(x = Value, fill = Scenario, color = Scenario)) +
geom_density(alpha = 0.5, size = 1) +
scale_x_continuous(labels = dollar) +
scale_fill_manual(values = c("Normal" = "#2c7fb8", "Alternative" = "#f03b20")) +
scale_color_manual(values = c("Normal" = "#2c7fb8", "Alternative" = "#f03b20")) +
labs(
title = "Monte Carlo Valuation Distribution Normal vs. Alternative Distributions",
subtitle = glue('Based on {comma(n_sim)} simulations.'),
x = "Firm Value",
y = "Density",
fill = "Scenario",
color = "Scenario"
) +
theme_minimal(base_size = 20)+
theme(legend.position = 'bottom')
Extension 2: Correlated Variables
In the baseline case, we have each random variable sampled independently from a normal distribution (i.e, i.i.d sampling). This can lead to higher variability and wider spread in the results that are unlikely to happen
By introducing a correlation structure between variables, we can constrain the possible outcomes, narrowing the distribution, leading to more realistic and stable simulations where changes in one variable influence the others
In what follows, we will keep \(r\) and \(g\) mean and standard deviation, but will now assume that the correlation between these two variables is \(\rho=0.9\):
- Whenever \(g\) is high, it might be because the firm is in its earlier stages, so \(r\) should be higher
- Whenever \(g\) is low, firm may have reached its market peak, so \(r\) should be lower
r_norm <- rnorm(n_sim, mean = 0.08, sd = 0.01)
g_norm <- rnorm(n_sim, mean = 0.03, sd = 0.005)
df_norm <- simulate_values(r_norm, g_norm, "Uncorrelated")
# Simulate r and g with negative correlation
mu_r <- 0.08
mu_g <- 0.03
sigma_r <- 0.01
sigma_g <- 0.005
corr_rg <- 0.9
# Covariance matrix
cov_matrix <- matrix(c(sigma_r^2, corr_rg * sigma_r * sigma_g, corr_rg * sigma_r * sigma_g, sigma_g^2),
nrow = 2, ncol = 2)
# Simulate r and g from the bivariate normal distribution
r_g_correlated <- mvrnorm(n_sim, mu = c(mu_r, mu_g), Sigma = cov_matrix)
r_corr <- r_g_correlated[, 1]
g_corr <- r_g_correlated[, 2]
df_corr <- simulate_values(r_corr, g_corr, "Correlated")
# Combine results
df_all <- rbind(df_norm, df_corr)
ggplot(df_all, aes(x = Value, fill = Scenario, color = Scenario)) +
geom_density(alpha = 0.6, size = 1.5) +
scale_x_continuous(labels = dollar, limits = c(0, 20000)) + # Set x-axis limits
scale_fill_manual(values = c("Uncorrelated" = "#2c7fb8", "Correlated" = "#f03b20")) +
scale_color_manual(values = c("Uncorrelated" = "#2c7fb8", "Correlated" = "#f03b20")) +
labs(
title = "Monte Carlo Valuation Distribution Uncorrelated vs. Strong Correlation",
x = "Firm Value",
y = "Density",
fill = "Scenario",
color = "Scenario"
) +
theme_minimal(base_size = 20)+
theme(legend.position = 'bottom')
Extension 3: Varying the Sample Size
In short, your Monte Carlo simulation outcome will depend on the number of simulations: the more simulations we run, the more stable and reliable the estimates become.
With smaller samples, (e.g., \(n<100\) or \(n<1,000\), results show high variability and may not capture the true distribution.
With larger samples, as \(n\) (e.g., \(100,000\) or \(1,000,000\)), the estimates converge toward the expected value
All in all, increasing sample size reduces random noise, improving the precision and the stability of simulation results
# Simulation sizes
n_sim_values <- c(100, 1000, 5000, 10000, 100000, 1000000)
# Store all results
all_results <- data.frame()
# Run simulations
for (n in n_sim_values) {
r_norm <- rnorm(n, mean = 0.08, sd = 0.01)
g_norm <- rnorm(n, mean = 0.03, sd = 0.005)
df <- simulate_values(r_norm, g_norm, comma(n))
all_results <- bind_rows(all_results, df)
}
all_results$Scenario <-factor(all_results$Scenario,levels=comma(n_sim_values))
# Plot histograms overlaid in one chart
ggplot(all_results, aes(x = Value, fill = Scenario)) +
geom_histogram(aes(y=..density..),bins = 60, color = "white") +
geom_density(size=0.5,fill=NA)+
scale_x_continuous(labels = scales::dollar,limits=c(0,25000)) +
facet_wrap(Scenario~.,scales='free')+
labs(
title = "Monte Carlo Valuation Distribution by Number of Simulations",
subtitle = "Varying the sample size for draws",
x = "Firm Value",
y = "Density",
fill = "Simulation Size"
) +
theme_minimal(base_size = 15) +
theme(legend.position = "bottom")
Extension 4: Bootstrapping
Bootstrapping is a resampling technique that draws samples with replacement from observed data
- It is useful to incorporate empirical variability without assuming a specific parametric distribution
- Helps capture the uncertainty present in historical or observed data
Given a sample \(x_1, x_2, \dots, x_N\), draw \(B\) bootstrap samples, you draw a given number of samples with replacement and compute the statistic of interest (e.g., mean, median) on each bootstrap sample
Bootstrapping can be used, for example, if we have a very limited series of historical values for our random values (e.g, \(n=20\)) but want to incorporate variability into the analysis
#Number of simulations
n_sim=10000
# Historical series
historical_cogs <- c(0.58, 0.6, 0.62, 0.59, 0.61,0.75,0.54,0.66,0.29,0.78)
historical_capex <- c(0.48, 0.5, 0.52, 0.49, 0.51,0.40,0.52,0.59,0.58,0.81)
# For each simulation, bootstrap 10 values and take the mean → length n_sim vector
boot_cogs_pct <- replicate(n_sim, mean(sample(historical_cogs, 10, replace = TRUE)))
boot_capex_pct <- replicate(n_sim, mean(sample(historical_capex, 10, replace = TRUE)))
simulate_values <- function(r, g, cogs_pct, capex_pct, label = 'Simulation') {
firm_value <- numeric(n_sim)
for (i in 1:n_sim) {
sales <- initial_sales
nwc_last <- sales * nwc_pct
fa_last <- sales * capex_pct[i]
fcf <- numeric(5)
for (t in 1:5) {
sales <- sales * (1 + g[i])
nwc_current <- sales * nwc_pct
fa_current <- sales * capex_pct[i]
delta_nwc <- nwc_current - nwc_last
delta_capex <- fa_current - fa_last
nwc_last <- nwc_current
fa_last <- fa_current
fcf[t] <- sales * (1 - cogs_pct[i]) - delta_capex - delta_nwc
}
fcf_terminal <- fcf[5] * (1 + g[i])
perp_value <- fcf_terminal / (r[i] - g[i])
discounted_fcf <- sum(fcf / (1 + r[i])^(1:5))
discounted_perp <- perp_value / (1 + r[i])^5
firm_value[i] <- discounted_fcf + discounted_perp
}
firm_value <- firm_value[is.finite(firm_value) & (r > g)]
data.frame(Value = firm_value, Scenario = label)
}
# Normal scenario (fixed percentages for COGS and CAPEX)
df_norm <- simulate_values(r_norm, g_norm, cogs_pct = rep(0.6,n_sim), capex_pct = rep(0.5,n_sim), label = "Normal")
# Bootstrap scenario (bootstrapped mean percentages)
df_bootstrap <- simulate_values(r_norm, g_norm, cogs_pct = boot_cogs_pct, capex_pct = boot_capex_pct, label = "Bootstrap")
#Plot
ggplot() +
geom_density(data = df_norm, aes(x = Value, fill = "Normal"), alpha = 0.5) +
geom_density(data = df_bootstrap, aes(x = Value, fill = "Bootstrap"), alpha = 0.5) +
scale_x_continuous(labels = dollar)+
labs(
title = "Monte Carlo Valuation Distribution: Normal vs. Bootstrap",
x = "Firm Value",
y = "Frequency",
fill = 'Approach',
color = "Scenario"
) +
theme_minimal(base_size = 20) +
theme(legend.position = "bottom")
References
Scheuch, Christoph, Stefan Voigt, and Patrick Weiss. 2023. Tidy Finance with R. Chapman & Hall/CRC. https://www.tidy-finance.org/r/.
Wickham, Hadley, Mine Cetinkaya-Rundel, and Garrett Grolemund. 2023. R for Data Science. O’Reilly Media. https://r4ds.had.co.nz/.